---
title: "7-5 Days alive and free of invasive or non-invasive ventilation"
description: |
This outcome was defined as the number of days alive and free of
invasive or non-invasive ventilation from randomisation to 28 days.
All patients dying within 28 days will be assigned zero free days.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
toc-depth: 5
freeze: true
---
```{r}
#| label: pkgs
#| code-summary: Load packages
library (ASCOTr)
library (tidyverse)
library (lubridate)
library (kableExtra)
library (patchwork)
library (cmdstanr)
library (posterior)
library (bayesplot)
library (ggdist)
library (lme4)
library (broom)
library (broom.mixed)
library (bayestestR)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
all_data <- readRDS (file.path (ASCOT_DATA, "all_data.rds" ))
# FAS-ITT
fas_itt_dat <- ASCOTr::: make_fas_itt_set (all_data)
fas_itt_nona_dat <- fas_itt_dat |>
filter (! is.na (out_dafv))
# ACS-ITT
acs_itt_dat <- ASCOTr::: make_acs_itt_set (all_data)
acs_itt_nona_dat <- acs_itt_dat |>
filter (! is.na (out_dafv))
# AVS-ITT
avs_itt_dat <- ASCOTr::: make_avs_itt_set (all_data)
avs_itt_nona_dat <- avs_itt_dat |>
filter (! is.na (out_dafv))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative" ), dir = "stan" )
ordmod <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_epoch_site" ), dir = "stan" )
ordmod_site <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_site" ), dir = "stan" )
logistic <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site_epoch" ), dir = "stan" )
```
```{r}
#| label: functions
#| code-summary: Functions
make_summary_table_anticoagulation <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (CAssignment = factor (
CAssignment,
levels = c ("C0" , "C1" , "C2" , "C3" , "C4" ),
labels = intervention_labels2 ()$ CAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_summary_table_antiviral <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (AAssignment = factor (
AAssignment,
levels = c ("A0" , "A1" , "A2" ),
labels = intervention_labels2 ()$ AAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (AAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Antiviral \n intervention ` = AAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
```
# Outcome Derivation
This outcome was defined as the number of days alive and free of positive pressure ventilation
(invasive or non-invasive) to 28 days. All patients dying within 28 days will be assigned zero free days.
The outcome is calculated for a patient as:
- missing, if day 28 mortality is unknown (`D28_PatientStatus` ) or if the patient was known to have been alive
at day 28 but number of days requiring ventilation (`D28_OutcomeDaysFreeOfVentilation` ) is unknown
- 0 if they died by day 28 (`D28_PatientStatus` )
- `min(28, D28_OutcomeDaysFreeOfVentilation)` otherwise
# Descriptive {#descriptive}
The overall distribution of days alive and free of ventilation (DAFV) are reported in @fig-dafv-dist for all participants in the ACS-ITT set.
There were few participants who required any ventilation during hospital who did not die (DAFV of 0),
```{r}
#| label: fig-dafv-dist
#| fig-cap: |
#| Distribution of days alive and free of ventilation
#| fig-subcap: FAS-ITT
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (dafv = ordered (out_dafv, levels = 0 : 28 ), .drop = F) %>%
mutate (p = n / sum (n))
p <- ggplot (pdat, aes (dafv, p)) +
geom_col () +
labs (
x = "Days alive and free of ventilation" ,
y = "Proportion"
)
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "outcome-7-5-dist-overall.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Anticoagulation
```{r}
#| label: tbl-7-5-summary-anticoagulation
#| code-summary: Summary of DAFH outcome by arm
#| tbl-cap: Summary of days alive and free of ventilation to day 28 by anticoagulation treatment group, FAS-ITT.
save_tex_table (
make_summary_table_anticoagulation (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-anticoagulation-summary" ))
make_summary_table_anticoagulation (fas_itt_dat)
```
```{r}
#| label: fig-dafh-dist-anticoagulation
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-anticoagulation.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-anticoagulation-avs-itt
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
pdat <- avs_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-anticoagulation-avs-itt.pdf" ), p, height = 3 , width = 6 )
p
```
## Antiviral
```{r}
#| label: tbl-7-5-summary-antiviral
#| code-summary: Summary of DAFH outcome by arm
#| tbl-cap: Summary of days alive and free of ventilation to day 28 by antiviral treatment group, FAS-ITT.
save_tex_table (
make_summary_table_antiviral (fas_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-antiviral-summary" ))
make_summary_table_antiviral (fas_itt_dat)
```
```{r}
#| label: fig-dafh-dist-antiviral
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-antiviral.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-antiviral-avs-itt
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention.
#| fig-height: 4
pdat <- avs_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = c ("Usual care" , "Nafamostat" )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.6 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
fpth <- file.path (pth, "outcome-7-5-descriptive-antiviral-avs-itt.pdf" )
ggsave (fpth, p + coord_flip (), height = 2.5 , width = 6 )
system (sprintf ("pdftoppm %s %s -png" , fpth, gsub (".pdf" , "" , fpth)))
p
```
```{r}
#| label: fig-dafh-dist-antiviral-acs-itt
#| fig-cap: |
#| Distribution of days alive and free of hospital by antiviral intervention.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (AAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Antiviral intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-antiviral-acs-itt.pdf" ), p, height = 3 , width = 6 )
p
```
## Age
```{r}
#| label: fig-dafv-dist-age
#| code-summary: DAFV by age
#| fig-cap: |
#| Distribution of days alive and free of ventilation by age group, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
dafv = fct_rev (ordered (out_dafv, levels = 0 : 28 )),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = dafv)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Sex
```{r}
#| label: fig-dafv-dist-sex
#| code-summary: DAFV by sex
#| fig-cap: |
#| Distribution of days alive and free of ventilation by sex.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Sex,
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (Sex) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Sex) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Sex, n)) +
geom_col () +
labs (
y = "Number of \n participants" )
p2 <- ggplot (pdat, aes (Sex, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Sex" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-sex.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Oxygen
```{r}
#| label: fig-dafv-dist-oxygen
#| code-summary: DAFV by oxygen
#| fig-cap: |
#| Distribution of days alive and free of ventilation by oxygen.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
supp_oxy2 = factor (supp_oxy2, labels = c ("Not required" , "Required" )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (supp_oxy2) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (supp_oxy2) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (supp_oxy2, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Supplemental oxygen" )
p2 <- ggplot (pdat, aes (supp_oxy2, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Supplemental oxygen" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-oxygen.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Country
```{r}
#| label: fig-dafv-dist-country
#| code-summary: DAFV by country
#| fig-cap: |
#| Distribution of days alive and free of ventilation by country.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = fct_infreq (Country),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (Country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (Country, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Country" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Site
```{r}
#| label: fig-dafv-dist-site
#| code-summary: DAFV by site
#| fig-cap: |
#| Distribution of days alive and free of ventilation by study site
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (Location),
dafv = ordered (out_dafv, levels = 0 : 28 )) %>%
complete (dafv = ordered (0 : 28 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Study site" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-dafv-dist-calendar
#| code-summary: DAFV by calendar date
#| fig-cap: |
#| Distribution of days alive and free of ventilation by calendar time.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Sample Cumulative Logits
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (CAssignment, out_dafv) %>%
complete (CAssignment, out_dafv, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafv) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafv != 28 )
ggplot (trt_logit, aes (CAssignment, rel_clogit)) +
facet_wrap ( ~ out_dafv) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
ggplot (trt_logit, aes (out_dafv, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (AAssignment, out_dafv) %>%
complete (AAssignment, out_dafv, fill = list (n = 0 )) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (AAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafv) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafv != 28 )
ggplot (trt_logit, aes (AAssignment, rel_clogit)) +
facet_wrap ( ~ out_dafv) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
ggplot (trt_logit, aes (out_dafv, rel_clogit)) +
facet_wrap ( ~ AAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
# Modelling
## FAS-ITT
```{r}
#| label: fit-model-fas-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
fas_itt_nona_dat,
ordmod,
outcome = "out_dafv" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save-fas-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects).
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-5-primary-model-fas-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-fas-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-5-primary-model-fas-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-5-primary-model-fas-itt-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
###### Diagnostics
```{r}
res$ fit$ summary (c ("alpha" , "beta" , "gamma_site" , "gamma_epoch" , "tau_site" , "tau_epoch" )) %>%
print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
###### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["alpha" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
### Posterior Predictive
```{r}
y_raw <- as.integer (levels (res$ dat$ y_raw))
y_ppc <- res$ drws$ y_ppc
y_ppc_raw <- rfun (\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols (fas_itt_nona_dat, tibble (y_ppc = y_ppc_raw))
grp_ppc <- function (grp, d = 0 : 27 ) {
ppc_dat %>%
group_by (grp = {{grp}}) %>%
summarise (
y_lteq = map_dbl (d, ~ mean (out_dafv <= .x)),
ypp_lteq = map (d, ~ rvar_mean (y_ppc <= .x)),
y_eq = map_dbl (d, ~ mean (out_dafv == .x)),
ypp_eq = map (d, ~ rvar_mean (y_ppc == .x))
) %>%
unnest (c (ypp_lteq, ypp_eq)) %>%
mutate (days = d) %>%
ungroup () %>%
pivot_longer (
y_lteq: ypp_eq,
names_to = c ("response" , "event" ),
names_sep = "_" ,
values_to = "posterior" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>%
filter (response == "ypp" , event == "eq" ),
aes (y = days)) +
facet_wrap ( ~ grp, nrow = 1 , scales = "free_x" ) +
stat_slabinterval (aes (xdist = posterior), fatten_point = 1 ) +
geom_point (data = dat %>% filter (response == "y" , event == "eq" ),
aes (y = days, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (
xlab, breaks = c (0 , 0.3 , 0.6 ),
sec.axis = sec_axis (
~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("Days" ) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_A <- grp_ppc (AAssignment)
pp_C <- grp_ppc (CAssignment)
pp_ctry <- grp_ppc (Country)
pp_epoch <- grp_ppc (epoch)
pp_site <- grp_ppc (site) %>%
left_join (ppc_dat %>% dplyr:: count (site, Country),
by = c ("grp" = "site" ))
p0 <- plot_grp_ppc (pp_A, "Antiviral" , "" )
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (
pp_site %>% filter (Country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (
pp_site %>% filter (Country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (
pp_site %>% filter (Country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (
pp_site %>% filter (Country == "NZ" ), "Sites New Zealand" , "" )
g1 <- (p0 | p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)
pth1 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-5-primary-model-fas-itt-ppc1.pdf" )
pth2 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-5-primary-model-fas-itt-ppc2.pdf" )
ggsave (pth1, g1, width = 6 , height = 5 , device = cairo_pdf)
ggsave (pth2, g2, width = 6 , height = 6 , device = cairo_pdf)
g1
g2
```
## ACS-ITT
```{r}
#| label: fit-model-acs-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
acs_itt_nona_dat,
ordmod,
outcome = "out_dafv" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save-acs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects).
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-5-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-5-primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-5-primary-model-acs-itt-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
## AVS-ITT
### Pre-specified
```{r}
#| label: fit-model-avs-itt-pre-spec
#| code-summary: Fit primary model
res <- fit_primary_model (
avs_itt_nona_dat |> mutate (ctry = droplevels (ctry)),
ordmod_site,
outcome = "out_dafv" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" , "ctry" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 1 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save-avs-itt-pre-spec
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-5-primary-model-avs-itt-summary-table-pre-spec" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-pre-spec
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-5-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_site_terms (
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-5-primary-model-site-terms-avs-itt-pre-spec.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Reduced Model
- excluding epoch, site, and country due to small sample size
```{r}
#| label: fit-model-avs-itt-reduced
#| code-summary: Fit primary model - AVS-ITT
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod0,
outcome = "out_dafv" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Required oxygen" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-5-primary-model-avs-itt-summary-table" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-5-primary-model-avs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
### Treatment Only
```{r}
res <- fit_primary_model (
avs_itt_nona_dat,
ordmod0,
outcome = "out_dafv" ,
vars = NULL ,
beta_sd_var = NULL ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
```
```{r}
#| label: odds-ratio-summary-table-avs-itt-trt-only
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table_rev (c (res$ drws$ AOR), "latex" ),
"outcomes/secondary/7-5-primary-model-avs-itt-summary-table-trt-only" )
odds_ratio_summary_table_rev (c (res$ drws$ AOR))
```
```{r}
#| label: fig-or-densities-avs-itt-model-trt-only
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-5-primary-model-avs-itt-odds-ratio-densities-trt-only.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```